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Abstract. We deduce the globular cluster formation history of the nearby elliptical galaxy, NGC 5128, by using a chemical 
enrichment model to accurately reproduce its observed metallicity distribution function (MDF). We derive the observed MDF 
using recently obtained U and B photometry of the NGC 5128 GC system, with (U — E) used as the metallicity indicator. 
Our results indicate that the GC system in this galaxy could be the product of two major GC formation episodes. The initial 
formation episode occured 1 1-12 Gyrs ago creating 65-75 percent of the mass in the GC system. This was followed by a second 
late formation episode which peaked 2-4 Gyrs ago and produced the remaining 25-35 percent of GC mass. 
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1. Introduction 

Globular cluster (GC) systems have been popular tools in de- 
ciphering the star formation histories (SFHs) of galaxies. Their 
ubiquitous presence in galaxies of all morphological types, 
coupled with evidence that their creation seems to accom- 
any major star formation episodes (e.g. lLarsen & RichtleJ 
1999h makes them usefu l tracers of galactic e v olution (e.g. 



Kissler-Patie et alJ Il998t Ivan den Berghl l2000t lYoon& Leej 
2002). In addition, GC systems are well represented by sim- 
ple stellar populations (SSPs) with stars of the same age and 
chemical composition, which makes them easy to study using 
SSP models (e.g. lYi et all2004 . 

The exact formation mechanism of GCs has been the 
subject of much recent debate. Many models for GC 
formatio n have been proposed including gase ous mer gers 
lAshman & ZepJ Il992h . in situ formation (e.g . lHarris et"al] 
1 1995b . multiphase collapse JForbes et alJll997l) , dissipation- 
less hierarchical merging (e.g. ICote et alJll998l EoOfl l2002h 
and hierarchical clustering ( Beasl ev eUuT200^i) /While none of 
these models can be conclusively excluded, the widespread dis- 
covery of multimodal metallicity distributions in GC popula- 
tions effectively rules out an extreme version of t he monol ithic 
collap se scenario for their formation (see e.g. iForbes et alJ 
Il997l) . Furthermore, the correlation between the mean metal- 
licity of GCs and galaxy luminosity indicates that chemical en- 
richment of the GC sy stem is intimately l inked to the evolu- 
tion of the host galax y ( Forbes et al. 1996; Durrell et al. 1996; 
IForbes & Fortd200ll1c6te et alJkOOOHjordan et alJl2004 . An 



elegant review of the GC format ion models mentioned above 
can be found in West et al. (2004). 

The object of this study is NGC 5128, the giant ellipti- 
cal galaxy in the nearby Centaurus Group (see llsraelH l998. 
for a comprehenive review of NGC 51 28). Located at a 
distance of approximately 3.6 Mpc (e.g. ISoria et alJ fl996>. 
NGC 5128 is the closest giant elliptical system with an es - 
timated GC population of 1550 + 350 (lHarris et alJ Il984h . 
It has been widely studied, not only because of its proxim- 
ity and relative brightness, but also because it displays un- 
usual physical features which suggest that this galaxy is a 
post-merger remnant. A pro minent dust lane containing young 
stars and HII regions (e.g. lUnger et al.| | 200(A IW ild & Eckart 
2000) and a se ries of optical shell s (Malin & Carted 1 198 3). 
which have H I ( Schiminovi ch et alJ ll994) and molecular CO 
JCharmandaris et alJl200 0) gas associated with them, are con- 
sidered strong evidence that NGC 5128 underwe nt a recent 
merge r event within the last 10 9 years. Recently. iReikuba et alJ 
(2004) suggested that star formation may have stopped as re- 
cently as 2 Myr ago in the north-eastern shell of NGC 5128. 
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formed a comprehensive study of the metallicity distribution 
of stars in the inner and outer halo of NGC 5128. In this 
study we deduce the formation history of the GC system of 
NGC 5128 by accurately reproducing its observed metallicity 
distribution function (MDF) using a chemical enrichment 
model. We derive the observed MDF of this elliptical galaxy 
using recently obtained U and B photometry of 210 clusters 
in its GC system JPeng et alJl2004al) . The integrated (U - B) 
colour is used as the metallicity indicator because it is sensitive 
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to metallicity via the opacity effect but relatively insensitive 
to the effective main sequence turn-off tempera ture (?Wf) an d 
therefore to age when T eff ~ 7000-12000 K JYi et aljboO^l) . 
Similar techniq ues u sing U band colour s have been used by 
iReikubal foOOli and Ijordan et alJ d2002t) who used (U - V) 
and Hubble Space Telescope WFPC2 (F336W - F547M), 
respectively. Although these colours are substantially better 
metallicity indicators than the previously used (V - /) or 
(B — V), they still change gradually with age and thus are not a s 
effective as (U - B) in determining metallicity jYi et al. 2004). 
Our work extends previous studies of NGC 5128 in a number 
of ways. The large number of confirmed GCs with U band 
photometry makes it possible to derive statistically significant 
results. In addition, the relative robustness of (U - B) as a 
metallicity indicator, compared to other optical colours, makes 
chemical enrichment an effective approach in modelling the 
formation history of the NGC 5128 GC system. The use of 
a consistent chemical enrichment code means that we can 
effectively transform metallicities into ages. 

We begin by briefly checking that, as might be expected 
from previous theoretical and observational results, a single 
starburst followed by passive evolution (a monolithic scenario) 
is incapable of reproducing the observed MDF of the NGC 
5128 GCs. Performing this check is not a redundant exercise 
because our modelling essentially yields relative likelihoods 
for various models to fit the observed MDF of NGC 5128. It 
is therefore instructive to compare the quality of fit between 
monolithic and non-monolithic scenarios. The main thrust of 
the paper, however, is to explore a double starburst scenario, 
analysing positions, timescales and relative strengths of the two 
star formation episodes that best explain the MDF of the NGC 
5128 GC system. We provide a coherent picture of the forma- 
tion history of NGC 5128 GCs based on double starburst sce- 
narios that give excellent fits to the observed MDF and show 
that our results are consistent with the spectroscopic study of 
IPeng et all J2004bl) . who use the age-sensitive Hp index to age- 
date the GCs in this dataset. 

2. Derivation of the observed NGC 5128 MDF 

We derive the observed MDF for the NGC 5 128 GC system by 
using photometric data in the U and B bands. We compute the 
metallicities by overlaying a model (U-B) vs. (B - V) grid on 
the GC data, with (U-B) used as the metallicity indicator. 

The metallicity of each GC is computed by identifying its 
two bounding iso-metallicity curves, followed by linear inter- 
polation between them. We do not use the same iso-metallicity 
curves to compute the metallicity error from the (U-B) er- 
rors. Instead we treat each error limit as a separate point on 
the grid and repeat the process described above. Given the 
irregular spacing and pseudo-horizontal nature of the model 
iso-metallicity curves (see Fig. 7 in Yi et al. 2004), this re- 
sults in a more accurate approximation of the metallicity errors 
for each GC. It should be noted that (U-B) colours (espe- 
cially redder colours) are affected by the age-metallicity de- 
generacy. The colour lying in a region of large degeneracy in- 
duces a larger metallicity error, since the corresponding metal- 
licity measurement becomes more uncertain. We also note that 
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Fig. 1. (U - B) colour distribution derived from final sample of 
GCs used in this study. 



lYi et alJ {2004) did not attempt to accurately derive the ages 
and metallicities of individual G Cs in thi s datase t. The age and 
metallicity derivations given in lYietalJ (2004) are rough es- 
timates derived from average values of (B - V) and (U-B) 
in broad bins. The method we use here provides far superior 
quantification of the metallicity by using a finely interpolated 
(U-B) grid. 

3. GC selection 

Our photometric data on NGC 5128 contains 210 GCs, some 
of which suffer from large uncertainties in the U and B magni- 
tudes. In this study we applied a cut of 0.1 mag to the uncer- 
tainty in (U - B), retaining 69 GCs out of the original sample 
of 210. The reasons for applying a rather stringent magnitude 
error cut are twofold. Firstly, some GCs have extremely large 
photometric errors. By using GCs with small error bars we are 
able to better quantify the metallicity of the GCs that we even- 
tually use along with their associated errors. This in turn results 
in more constrained estimates for the characteristics of the star 
formation episodes. Secondly, having performed this magni- 
tude cut we find that the MDF for the entire sample of 210 GCs 
is very similar to o ur reduced sample; a Kolmogorov-Smirnov 
(KS) test (see e.g. Wall 1996) between the entire sample and 
the distribution derived after the uncertainty cut yields a KS 
probability of 42 percent, indicating that the two samples are 
consistent at the 95 percent confidence level. 15 percent of the 
final sample of 69 GCs lie within a projected radius of 6' from 
the centre of the galaxy (NGC 5128 has an apparent size of 18' 
x 14') compared to 25 percent of the original sample. 56 per- 
cent of the final sample lies within 10' compared to 59 percent 
of the original sample. The faint end of the GCLF of the final 
sample lies at V — 19.7 mag, while the faint-magnitude cuts 
are 20.6 and 20.4 in the U and B bands respectively. 

Our aim is to maximise the accuracy of our analysis. The 
magnitude cut we employ increases the robustness of our final 
solutions, without transforming the original distribution signif- 
icantly. Fig.^presents the (U - B) colour distribution derived 
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Fig. 2. TOP PANEL: Metallicity distribution derived from ( U - 
B) colours of 69 GCs which have (U - B) photometric errors 
less than 0.1 mag. All photometry is corrected for Galactic ex- 
tinction. BOTTOM PANEL: Comparison between the derived 
metallicity distribution (shown in top panel) and that derived 
by Harris et al. (2003) using the metallicity sensitive (C - T\) 
index. 



from our final sample of clusters and Fig.|2presents the metal- 
licity distribution derived from the (U — B) colours of our final 
sample of clusters. We note that our derived MDF is consistent 
with the NGC 5128 MDF derived using the (C - T{) index for 
205 clusters in lHarris et alJ J2003h . This comparison is shown 
in the bottom panel of Fig. [2 

4. Star formation and gas infall prescriptions 

Feedback processes from star formation remain poorly un- 
derstood. Supernovae typically inject large amounts of en- 
ergy into the inter-stellar medium. The energy deposited can 
be either kinetic or thermal and the exact description of 
thes e processes remains the subject of much debate (see 
e.g. lKadll992t iMihos & Hernauisdfl994l lYeoes etai]ll997l 
iFerreras et alJl2002t ISilkll2003h . Such feedback processes can 
interrupt or accelerate star formation, so the simple coupling 
between star formation and the gas density implied by, for ex- 
ample, a Schmidt law (Schmidt 1959) may not hold. In this 
model, we decouple the star formation rate from the gas den- 
sity in the system to qualitatively account for such feedback 
processes from star formation. We take the star formation rate 
to be a sequence of gaussian functions independent of the gas 
density 



¥0) = 2 



fn 



V27T0-,, 



exp ■ 



~(t - T„) 2 



2<rl 



(1) 



In Eq.^/, is the normalisation amplitude, <x„ the temporal 
spread and t„ the centre-point of starburst n. With the gaussian 
starbursts normalised as shown above, we keep the total area 
under T'(f) and hence the total stellar mass equal to 1 by ensur- 
ing that £ii f n = 1. We checked that our results do not depend 
strongly on our choice of mathematical function by repeating 



our analysis using lorentzian starbursts, and find that our con- 
clusions remain unchanged. 

We assume that the system begins with no gas and use an 
exponentially decaying gas infall rate 



gin(t) = — exp(— ). 

'in V * in ' 



(2) 



Here r,„ is the infall timescale, and 77 is a parameter that con- 
trols the final gas fraction in the system and is defined below. 
Noting that, 
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where ty is the age of the universe in Gyrs, we can force the 
system to have a final gas fraction gf rac by choosing 



1 



gfrac 



so that, 



n = 



1 



■gfr, 



(4) 



(5) 



We fix the infall timescale r,„ to 5 Gyrs, which leaves the fol- 
lowing 3n free parameters in the model for an n starburst sce- 
nario 

• ■ • ,fn-2,T\, ■ ■ ■ ,T„-U CT 1, . . . ,(T n -i, gfrac). 

We note that due to the decoupling of the star formation rate 
from gas density, changing the value of the r,„ has no impact 
on our subsequent analysis. Essentially, the behaviour of the 
system remains unchanged as long as there is a reservoir of gas 
available to make stars. 

Given the prescriptions given above for the star formation 
and gas infall rates, we compute the metallicity evolution of the 
GCs following the s tandard set of chemica l enrichment equa- 
ti ons as describ ed in Ferreras & Silk (2000). We use the yields 
of Thielemann et al. ( 1996) for stellar masses M* > 1OM and 
Ivan den H oek & Groenewegenl i 1 997h for lower mass stars. For 
a given set of parameters we integrate the equations to allow us 
to compute the MDF, that is then compared with the observed 
MDF via a KS test. 

5. Derivation of best-fit parameters 

We begin our analysis by comparing a single starburst sce- 
nario, which describes a monolithically evolving system, to 
the more realistic ones that underpin this study. We then re- 
fine our analysis by adding another starburst and perform a de- 
tailed study of the resulting six dimensional parameter space. 
We also check that small perturbations to the two major star 
formation episodes do not affect the best-fit values for our pa- 
rameters. This sequence is repeated for four final gas fractions - 
gfrac = 0.1,0.2,0.35,0.5. Since we are modelling an elliptical 
galaxy, the value of gf rac is likely to be low, probably some- 
where in the region 0. 1 -0.35. In a recent work. lSanderson et alJ 
(2003) finds the (hot) gas fraction in ellipticals within R200 to 
be approximately 0. 1 . 
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We compare the MDF predicted by our model to the ob- 
served GC MDF using the KS test. The best-fit values of pa- 
rameters are calculated by marginalising their KS probabili- 
ties, which involves summing (integrating) out the probabil- 
ity dependence of all other parameters, leaving a probability 
distribution for each parameter that is independent of all oth- 
ers in the model. For example, in a double starburst scenario, 
the marginalised probability distribution for the parameter n is 
calculated as 



T2i, CT\ /, CT2k), 



(6) 



where g is a normalising constant which ensures that P(t-[) in- 
tegrates to 1 . This probability distribution will not, in general, 
be symmetric and therefore should not be approximated by a 
gaussian. 

We take the best-fit value of a parameter as the value at 
which its marginalised probability function peaks. In addition, 
we define errors for these parameters superficially similar to 
the one-sigma errors frequently used for gaussian probability 
distributions. If X is the most likely value of a parameter x, 
then we define the positive error x + as 



P(x)dx = 0.34 I P(x)dx, 
x Jx 

and similarly the negative error, x_, as 
I P(x)dx = 034 I P(x)dx. 

%J X- %J—(X 



(7) 



(8) 



Finally, given the fact that each GC has a metallicity error 
associated with it, we check that the positions of the most likely 
values for our model parameters are robust with respect to these 
errors. We perform this check by adding random noise to our 
observed MDF and producing KS values for these noisy MDFs. 
The noise is added by running through each GC and changing 
its metallicity to one of its error points randomly. This process 
is repeated a hundred times to generate a hundred noisy MDFs. 
We have checked that the KS values do not deviate significantly 
from the fiducial value and find that the KS values are suitably 
robust (within 6 percent of the original values) when the fidu- 
cial MDF is perturbed in this way. This is mainly a result of 
choosing GCs with small photometric errors. 

6. Results and analysis 

6.1. Double starburst scenarios 

For a single starburst scenario, the maximum KS probability we 
achieve in our model is around 2 percent, significantly lower 
than the double starburst scenarios we explore later. The most 
favoured values of r and cr for the single starburst scenario are 
1 Gyr and 1 .35 respectively. A single starburst model fails to 
reproduce, with any acceptable accuracy, the observed MDF of 
NGC 5128. 

With two starbursts we obtain a six dimensional parameter 
space with the star formation rate as the sum of two gaussians 



fx -(t-rQ 2 , -(f-r 2 ) 2 

exp — — — + — exp — — — . (9) 



V27TO-, 



2crf 



V27TCT 2 



2c\ 



Table 1 shows the best-fit values for the starburst parameters 
(fx, fi, fx> i"2, o"ti 0"2) f° r final gas fractions of 0.5 to 0.1. Note 
that in Table 1, Subscript '1' indicates the initial (early) star- 
burst and Subscript 2 indicates the second (late) starburst. 

In Fig. [3] we present the best KS values produced in the 
double starburst models. The single starburst case is shown, 
for comparison, by the solid line at 2 percent. In Fig. 0] we 
present three GC formation histories which yield the highest 
KS probabilities (see shaded region in Fig. [3} by showing the 
formation histories for favoured values (/! , f-i) in a system with 
a final gas fraction of 0.35. 

Fig. [3] suggests that an intermediate gas fraction of 0.35 
with (fx,fi) - (0.35,0.65) is the most favoured model, so that 
roughly one-third of the mass in the GC system is predicted 
to have been formed in a late starburst centred around 2 Gyrs 
ago. Fig. [5] shows the mass fractions formed at various ages 
and metallicities in this model, and Fig.[6]compares the MDF 
predicted by this model to the observed MDF. 

Although the favoured model provides the best-fit, the 
shaded region in Fig. [3] shows a set of degenerate solutions 
consisting of models producing high KS values, all of which 
reasonably reproduce the observed NGC 5 128 MDF and there- 
fore cannot be completely discounted. Changing the relative 
sizes of the starbursts even further results in the KS probabil- 
ity declining rapidly. This is clearly expected, since in the limit 
f\ — > 0, we recover a single starburst, already established as 
being incapable of reproducing the observed NGC 5128 MDF. 
While the position of the first starburst is fairly robust in the 
best-fit models, the second starburst tends to vary, so that we 
are left with a dispersion of ~ 2 Gyrs in the second starburst's 
position. This reflects the fact that uncertainties in colours and 
in the age-metallicity degeneracy (discussed in Sect. 2) induce 
uncertainties in the metallicity derivations for our GCs, which 
in turn limits our resolution of the position of the late GC for- 
mation episode. Averaging over all the best-fit models (shown 
in the shaded region in Fig. [3J and taking the one-sigma er- 
rors into account for the parameters given in Table 1, our re- 
sults therefore indicate that the observed MDF of the NGC 
5128 GCs can be reproduced accurately by a formation sce- 
nario where roughly 65-75 percent of the mass in the GC sys- 
tem was formed in an initial starburst 11-12 Gyrs ago, followed 
by a second late starburst which peaked 2-4 Gyrs ago and pro- 
duced the remaining 25-35 percent of GC mass. 

6.2. Perturbations to the double starburst case 

As mentioned above, the n starburst scenario contains 3n free 
parameters. Clearly, as we make more free parameters avail- 
able to the model, we should expect to do better in terms of 
matching the observed NGC 5128 MDF. We note first that the 
KS values from the double starburst scenarios are significantly 
large. Refining the model further would not make an apprecia- 
ble difference in the KS probabilities, certainly not as much as 
in the transition from the single to the double starburst case. In 
addition, given the errors in the observed MDF and the resolu- 
tion of our model, further refinement would not provide much 
greater insight into the GC formation history of this galaxy. 
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Table 1. Best-fit values and errors for the double starburst scenarios with final gas fractions in the 
range 0.1 - 0.5. Note that Subscript '1' indicates the initial (early) starburst and Subscript 2 the 
second (late) one. 
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Fig. 3. Best KS probabilities for various values of fi and for four fi- 
nal gas fractions g/ rac . The shaded region shows a set of degenerate 
models that have high KS probabilities and therefore cannot be dis- 
counted completely. The KS probability for a single starburst scenario 
is shown by the solid line at 2 percent. 



Fig. 5. Age-metallicity grid showing the fraction of stellar mass 
predicted by the favoured model as a function of age and metal- 
licity. The key indicates stellar fractions corresponding to the 
levels of shading in the plot. 
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Fig. 4. Three GC formation histories corresponding to the best KS 
probabilities shown in the shaded region in Fig. [5] 
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Fig. 6. Comparison between observed MDF and that predicted 
by the favoured model (which gives the highest KS probabil- 

lity). 



Therefore, instead of repeating the study above with a 
wholesale three starburst analysis, we chose to treat the third 
starburst as a perturbation on the double starburst system. We 
were interested in studying how the KS probabilities change if 
we apply a small perturbation to each of the major star forma- 
tion episodes. This simple perturbation analysis gave an indica- 
tion of how localised these double starburst KS hotspots were 
in the nine parameter triple starburst space. 

Fig. shows the best KS values achieved after applying a 
5 percent pertubation i.e. f perturb - 0.05 as a fraction of the 
best KS values in the corresponding unperturbed double star- 
burst scenario. In almost all cases the KS value improves only 
marginally. We conclude from this analysis that the KS hotspot 
is not extremely localised in the nine parameter triple starburst 
space and does appear to be stable at least under perturbations 
of sizes 5 percent and below. Although we have sampled only 
a small fraction of this nine-dimensional parameter space, the 



resistance to small perturbations indicates that the double burst 
scenarios are reasonably stable and robust solutions. 

7. Discussion and conclusions 

We studied the formation history of the NGC 5128 GC system 
by using a chemical enrichment model to accurately reproduce 
its observed MDF derived from recently obtained U and B pho- 
tometry. Our results show that the GC system in this galaxy 
could be the product of two major GC formation episodes, al- 
though we do not have adequate resolution in the model or in 
the observed MDF to make any reliable claim about the possi- 
bility of additional formation events. 

The double starburst analysis produces high KS probabil- 
ities and therefore good fits to the observed MDF, with best- 
fit models (shaded region in Fig. [5J that favour an initial GC 
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Perturbation size= 5 percen 



g frao =0.20 

g Frac =0.35 

g Frnc =0.50 



Fig. 7. KS values for perturbed double starburst scenarios as a fraction 
of the KS values in the unperturbed case. The size of the perturbation 
applied to the system is 5 percent. 



formation episode 1 1-12 Gyrs ago which produced 65-75 per- 
cent of mass in the GC system, with a second late formation 
episode approximately 2-4 Gyrs ago producing the remaining 
25-35 percent. The late starburst results in a small fraction (~ 5 
percent) of the stellar mass in the GC system potentially having 
ages less than 1 Gyr. If NGC 5128 has over 1500 clusters, then 
we might expect a handful of those clusters to have ages less 
than 1 Gyr. We have not made any a priori assumptions about 
the driving mechanisms behind star formation episodes in our 
model. While a single starburst followed by passive evolution 
can clearly be discounted, the nature of our model admits any 
scenario where multiple episodes of GC formatio n are possible. 
Our results are in general agreement with iBeaslev et alJ 



(2003), who obtain good fits to the observed NGC 5128 MDF 
using a semi-analytical galaxy formation model with metal- 
poor GC formation halted at z ~ 5. A quantitative similarity 
between this study and that of Beasle v"et alJ J2003I) is that the 
initial GC formation episode (first gaussian starburst) in our 
model, which gave rise to the metal-poor GC subpopulation, 
decays rapidly to virtually zero (see Fig.|4j within ~ 2 Gyrs of 
the formation of the galaxy; i.e. there is effective truncation of 
metal-poor GC production at high redshift (z ~ 3 — 4). This 
closely resembles the truncation employed by IBeaslev et alJ 
(2003), and the similarity is probably due to the fact that our 
model is comparable to the chemical en richment prescriptio n 
employed in the semi-analytical model of lBeaslev et al.l (12003). 
Since the object of both studies is to reproduce metallicity dis- 
tributions this resemblance is not unexpected. 

This study would not be comple te witho ut a comparison 
of our results to those of IPeng et al.l ll2004bll . who performed 
a spectroscopic analysis of the GC dataset used in this pa- 
per. Based on the age-sensitive Hp index, they conclude that 
metal-poor GCs in NGC 5 128 have ages com parable t o those in 
the Milky Way, i.e. ~ 12 Gyrs (e.g lKrauss & Chaboverl2003l) . 
and that the metal-rich GCs are consistent with a mean age of 
5^ Gyrs. We thus find that our age estimates, derived using a 
chemical enrichment approach to exploit metallicity-sensitive 



photometric colours, are consistent with a study of the same 
objects using age-sensitive spectroscopic indices. 

Our study demonstrates the potential of a chemical enrich- 
ment approach in deciphering the formation histories of the 
GC system in galaxies. As more age and metallicity-sensitive 
spectro-photometric data become available, methods such as 
the one used in this study will enable us to set robust constraints 
on the way galaxies either form or incorporate GCs, crucial to 
our understanding of galaxy formation. 
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